load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

;location_f="/your_path_here"
location_f="/Users/ivana/WORK_BSC/data_to_publish/CCSM4som/"

filsw= systemfunc ("ls "+location_f+"lowice*2dvars*.nc")
filsc= systemfunc ("ls "+location_f+"ctrl*2dvars*.nc")

fw    = addfiles (filsw, "r")
fc    = addfiles (filsc, "r")

ListSetType (fw, "join")
ListSetType (fc, "join")
print(fw)
print(fc)

;*******************************************;
;Read in variables
;*******************************************
lat=fc[0]->lat
lon=fc[0]->lon
time=fc[0] ->time
rad  = 4.0 * atan(1.0) / 180.
clat=dble2flt(cos(lat*rad))

nlat  = dimsizes(lat)
nlon  = dimsizes(lon)
ntime= dimsizes(time)
print(ntime)

lhflxc=fc[:]->LHFLX
lhflxw=fw[:]->LHFLX

shflxc=fc[:]->SHFLX
shflxw=fw[:]->SHFLX

fsnsc=fc[:]->FSNS
fsnsw=fw[:]->FSNS

flnsc=fc[:]->FLNS
flnsw=fw[:]->FLNS

fsntc=fc[:]->FSNT
fsntw=fw[:]->FSNT

flntc=fc[:]->FLNT
flntw=fw[:]->FLNT

tempsc=flntc
tempsc=fsntc-flntc+( -fsnsc+flnsc+lhflxc+shflxc)

tempsw=flntw
tempsw=fsntw-flntw +(-fsnsw+flnsw+lhflxw+shflxw)

printVarSummary(tempsc)
printVarSummary(tempsw)


;************************************************
; average over the simulations/years
;************************************************

tcs=tempsc(0,:,:,:)
tcs=dim_avg(tempsc(time|:,lat|:,lon|:,ncl_join|:))
tws=tempsw(0,:,:,:)
tws=dim_avg(tempsw(time|:,lat|:,lon|:,ncl_join|:))

finfo_tcs=tcs(0,:,:)
finfo_tcs=dim_avg(tcs(lat|:,lon|:,time|:))
finfo_tws=tws(0,:,:)
finfo_tws=dim_avg(tws(lat|:,lon|:,time|:))


finfo_tcs_z=finfo_tcs(:,0)
finfo_tcs_z=dim_avg(finfo_tcs(lat|:,lon|:))
finfo_tws_z=finfo_tws(:,0)
finfo_tws_z=dim_avg(finfo_tws(lat|:,lon|:))

;************************
; calculate the anomaly:
;************************

  finfo_ano=finfo_tws_z
  finfo_ano=finfo_tws_z- finfo_tcs_z
  finfo_ano=finfo_ano*clat

  finfo_ano!0="lat"
  finfo_ano&lat=lat
  print(finfo_ano)

;________________
;write data

;name ="columnnet_ano_zonal_ccsm4som"
;outf = addfile(name+".nc","c")
;outf->time = time
;outf->lat  = lat
;outf->finfo_ano_z   = finfo_ano

;****************************************
; Create plot
;****************************************
wks = gsn_open_wks("pdf","zonal_columnnetflux_ccsm4som")

   plot    = new (1,"graphic")

   res                      = True          ; individual plot
   res@gsnDraw              = False
   res@gsnFrame             = False
   res@xyLineThicknessF     = 2.0
   res@xyMonoLineColor         =True
   res@xyLineColor         = "black"
   res@gsnYRefLine           = 0.0             ; create a reference line
   res@gsnAboveYRefLineColor = "red"              ; above ref line fill red
   res@gsnBelowYRefLineColor = "blue"             ; below ref line fill blue

   res@xyDashPattern        = 0                  ; Make curves all solid

   res@vpHeightF= 0.4                    ; change aspect ratio of plot
   res@vpWidthF = 0.8

   res@trYMinF=-1
   res@trYMaxF=1.5

   res@tiXAxisFontHeightF = 0.010
   res@tiYAxisFontHeightF = 0.015

   res@gsnLeftString = ""
   res@tiYAxisString = "Column net flux ano"
   res@tiXAxisString = ""

   plot(0)  = gsn_csm_xy (wks,lat,finfo_ano,res)
   res@tiYAxisString = "Column net flux ano"

;***** make panel *****


   resP                     = True          ; panel resources
   resP@gsnMaximize         = True

   resP@txFontHeightF = 0.015
   gsn_panel(wks,plot,(/1,1/),resP)

